In [1]:

    

%load_ext snakeviz
%autoreload 2
import importlib


    

In [2]:

    

import sys
sys.path.append('/Users/dsuess/Code/Pythonlibs')

import mpnum as mp
from mptomo.tiht import tiht_np, tiht_mpa
import mptomo.tiht as tiht
import numpy as np
import h5py

from numpy import float_
from tools.helpers import Progress, Timer, TimelyProgress, RuntimeSlice, CountSlice
from itertools import islice


    

In [10]:

    

RGEN = np.random.RandomState(seed=0)
SITES = 6
RANK = 4
LDIM = 5
MEAS = 6 * SITES * RANK**2 * LDIM
BATCHSIZE = 150
RETINFO = ['mu']

X = mp.random_mpa(SITES, LDIM, RANK, randstate=RGEN,
                  normalized=True, dtype=float_)
X /= mp.norm(X)

A = [mp.random_mpa(SITES, LDIM, 1, randstate=RGEN, normalized=False, dtype=np.float_) / MEAS
     for _ in range(MEAS)] 

y = np.array([mp.inner(a, X) for a in A])


    

In [11]:

    

stepfun = tiht.memory_stepsize(tiht.locally_projected(tiht.adaptive_stepsize()),
                               const_steps=50)
solution = tiht_mpa(A, y, 2 * RANK, batchsize=BATCHSIZE, stepsize=stepfun, retinfo=['mu'])
result = list(a for a in Progress(RuntimeSlice(solution, 60), rettime=True))


    

    




---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-11-b29c6a0e9ebb> in <module>()
      2                                const_steps=50)
      3 solution = tiht_mpa(A, y, 2 * RANK, batchsize=BATCHSIZE, stepsize=stepfun, retinfo=['mu'])
----> 4 result = list(a for a in Progress(RuntimeSlice(solution, 60), rettime=True))

<ipython-input-11-b29c6a0e9ebb> in <genexpr>(.0)
      2                                const_steps=50)
      3 solution = tiht_mpa(A, y, 2 * RANK, batchsize=BATCHSIZE, stepsize=stepfun, retinfo=['mu'])
----> 4 result = list(a for a in Progress(RuntimeSlice(solution, 60), rettime=True))

/Users/dsuess/Code/Pythonlibs/tools/helpers.py in __iter__(self)
    134         """Fetch next object from the iterable"""
    135         self.start()
--> 136         for runtime, val in self._iterable:
    137             self.update(min(runtime, self._iterable.runtime))
    138             yield val if not self._rettime else (runtime, val)

/Users/dsuess/Code/Pythonlibs/tools/helpers.py in __iter__(self)
     94     def __iter__(self):
     95         starttime = time.time()
---> 96         for val in self._iterable:
     97             runtime = time.time() - starttime
     98             yield runtime, val

/Users/dsuess/Code/MPO Tomography/mptomo/tiht.py in tiht_mpa(A, observation, bdim, X_init, batchsize, retinfo, stepsize, compargs)
    219         g.compress(bdim=2 * bdim, **compargs)
    220 
--> 221         mu = stepsize(A, g, X_for_stepsize)
    222 
    223         X += float(mu) * g

/Users/dsuess/Code/MPO Tomography/mptomo/tiht.py in stepsize(A, g, X)
    159     def stepsize(A, g, X):
    160         if stepsize.counter >= const_steps:
--> 161             mu = stepfun(A, g, X)
    162             stepsize.memory.append(mu)
    163             stepsize.counter = 1

/Users/dsuess/Code/MPO Tomography/mptomo/tiht.py in stepsize(A, g, X)
    148     def stepsize(A, g, X):
    149         # its actually not a projection, but it works surprisingly
--> 150         proj = mp.MPArray([_local_contraction(l)[None, ..., None] for l in X[:-1]])
    151         g = mp.partialdot(proj, g, start_at=0)
    152         g.compress(**compargs)

TypeError: 'MPArray' object is not subscriptable

In [ ]:

    

stepfun = tiht.memory_stepsize(tiht.fully_projected(tiht.adaptive_stepsize()), 
                               const_steps=50)
solution = tiht_mpa(A, y, 2 * RANK, batchsize=50, stepsize=stepfun, retinfo=['mu', 'X_for_stepsize'])
result2 = list(a for a in Progress(RuntimeSlice(solution, 60), rettime=True))


    

In [113]:

    

pl.plot([info['mu'] for _, (_, info) in result])
pl.plot([info['mu'] for _, (_, info) in result2])


    

    
Out[113]:





[<matplotlib.lines.Line2D at 0x10f719cf8>]






    

In [114]:

    

pl.plot([t for t, _ in result], [mp.normdist(X, Xsharp) for _, (Xsharp, _) in result], label='1')
pl.plot([t for t, _ in result2], [mp.normdist(X, Xsharp) for _, (Xsharp, _) in result2], label='2')

pl.legend()
pl.show()

pl.plot([mp.normdist(X, Xsharp) for _, (Xsharp, _) in result], label='1')
pl.plot([mp.normdist(X, Xsharp) for _, (Xsharp, _) in result2], label='2')

pl.legend()
pl.show()


    

In [64]:

    

Xs = [info['X_for_stepsize'] for _, (_, info) in result2]


    

In [66]:

    

X = Xs[5]
X.normalize(left=len(X) - 1)

# get the left-eigenvectors
Us, _ = X.split(len(X) - 2)
Us = Us.reshape([s + (1,) for s in Us.pdims[:-1]] + [Us.pdims[-1]])
proj = mp.dot(Us.conj(), Us, axes=(1, 1))


    

In [74]:

    

a = np.array([1, 2, 3])


    

In [82]:

    

def local_proj(l):
    l = np.rollaxis(l, 1).reshape((l.shape[1], -1))
    return np.tensordot(l.conj(), l, axes=(1, 1))


    

In [95]:

    

def locally_projected(stepfun, compargs={}):
    def stepsize(A, g, X):
        proj = mp.MPArray([local_proj(l)[None, ..., None] for l in X[:-1]])
        g = mp.partialdot(proj, g, start_at=0)
        g.compress(**compargs)

        return stepfun(A, g, _)
    return stepsize


    

In [90]:

    

[local_proj(l)[None, ... , None].shape for l in X[:-1]]


    

    
Out[90]:





[(1, 10, 10, 1), (1, 10, 10, 1)]

In [91]:

    

proj = mp.MPArray([local_proj(l)[None, ..., None] for l in X[:-1]])


    

In [92]:

    

mp.normdist(mp.dot(proj, proj), proj)


    

    
Out[92]:





6.0975218522976471

In [48]:

    

def overlap(mpa1, mpa2):
    return mp.inner(mpa1, mpa2) / mp.norm(mpa1.copy()) / mp.norm(mpa2.copy())

def get_overlaps(result, *args, **kwargs):
    for _, (_, info) in result:
        g = mp.sumup(info['A_'], weights=(info['y_'] - info['y_sharp']))
        compr, _ = g.compression(*args, **kwargs)
        yield overlap(g, compr), compr.bdim


    

In [49]:

    

axl = pl.gca()
axr = axl.twinx()

with_relerr = list(get_overlaps(Progress(result), 'svd', relerr=0.3))
axl.plot([ol for ol, _ in with_relerr])
axr.plot([bdim for _, bdim in with_relerr], ls=":")

with_relerr = list(get_overlaps(Progress(result), 'svd', bdim=4 * RANK))
axl.plot([ol for ol, _ in with_relerr])
axr.plot([bdim for _, bdim in with_relerr], ls=":") 

with_relerr = list(get_overlaps(Progress(result), 'var', bdim=4 * RANK))
axl.plot([ol for ol, _ in with_relerr])
axr.plot([bdim for _, bdim in with_relerr], ls=":")


    

    




 98% (61 of 62) |######################## | Elapsed Time: 0:00:25 ETA:  0:00:00





    
Out[49]:





[<matplotlib.lines.Line2D at 0x11140bda0>]






    

In [29]:

    

np.random.shuffle?


    

In [30]:

    

result = dict()

for batchsize in reversed([50, 100, 150, 200, 300, 400, 500, 600, 1000]):
    solution = tiht_mpa(A, y, 2 * RANK, batchsize=batchsize,
                       stepsize={'method': 'adaptive_noproj', 'recalc': 50, 'join': 'average', 'scale': 1.3})
    result[batchsize] = list(a for a in Progress(RuntimeSlice(solution, 30), rettime=True))


    

    




100% |#######################################| Elapsed Time: 0:00:30 / 00:00:30
100% |#######################################| Elapsed Time: 0:00:30 / 00:00:30
100% |#######################################| Elapsed Time: 0:00:30 / 00:00:30
100% |#######################################| Elapsed Time: 0:00:30 / 00:00:30
100% |#######################################| Elapsed Time: 0:00:30 / 00:00:30
100% |#######################################| Elapsed Time: 0:00:30 / 00:00:30
100% |#######################################| Elapsed Time: 0:00:30 / 00:00:30
100% |#######################################| Elapsed Time: 0:00:30 / 00:00:30
100% |#######################################| Elapsed Time: 0:00:30 / 00:00:30

In [24]:

    

with h5py.File('batchsize_scan_8.h5', 'w') as buf:
    for varname in ['SITES', 'RANK', 'LDIM', 'MEAS']:
        buf.attrs[varname] = globals()[varname]
    
    X.dump(buf.create_group('X'))
    X.dump(buf.create_group('A'))
    
    for key, val in result.items():
        group = buf.create_group(str(key))
        for time, X_sharp in val:
            X_sharp.dump(group.create_group(str(time)))


    

In [25]:

    

for key, val in result.items():
    pl.plot([t for t, _ in val], [mp.normdist(X, Xsharp) for _, Xsharp in val], label=key)

pl.legend()
pl.show()

for key, val in result.items():
    pl.plot([mp.normdist(X, Xsharp) for _, Xsharp in val], label=key)

pl.legend()
pl.show()


    

In [26]:

    

del result


    

In [12]:

    

solution = tiht_mpa(A, y, 2 * RANK, batchsize=150, 
                    retinfo=['mu'])
result = list(a for a in Progress(RuntimeSlice(solution, 30), rettime=True))
pl.plot([mp.normdist(X, X_sharp) for _, (X_sharp, _) in result])
pl.gca().twinx().plot([info['mu'] for _, (_, info) in result], color='b')


    

    




100% |#######################################| Elapsed Time: 0:00:30 / 00:00:30






    
Out[12]:





[<matplotlib.lines.Line2D at 0x10fab6048>]






    

In [20]:

    

%%snakeviz

solution = tiht_mpa(A, y, 2 * RANK, batchsize=150, retinfo=['mu'],
                    stepsize={'method': 'adaptive_noproj', 'recalc': 50, 'join': 'average', 'scale': 1.3})
result_r = list(a for a in Progress(RuntimeSlice(solution, 30), rettime=True))


    

    




100% |#######################################| Elapsed Time: 0:00:30 / 00:00:30






    




 
*** Profile stats marshalled to file '/var/folders/25/5c0947xn1md22_9d8vgnm45w0000gn/T/tmpwwvjqsng'. 

In [26]:

    

solution = tiht_mpa(A, y, 2 * RANK, batchsize=150, retinfo=['mu'],
                    stepsize={'method': 'adaptive_noproj', 'recalc': 50, 'join': 'average', 'scale': 1.3},
                    compargs={'method': 'svd'})
result_r = list(a for a in Progress(RuntimeSlice(solution, 30), rettime=True))

pl.plot([mp.normdist(X, X_sharp) for _, (X_sharp, _) in result])
pl.plot([mp.normdist(X, X_sharp) for _, (X_sharp, _) in result_r])

ax2 = pl.gca().twinx()
ax2.plot([info['mu'] for _, (_, info) in result], ls="--")
ax2.plot([info['mu'] for _, (_, info) in result_r], ls="--")


    

    




100% |#######################################| Elapsed Time: 0:00:30 / 00:00:30






    
Out[26]:





[<matplotlib.lines.Line2D at 0x113db2b70>]






    

In [10]:

    

result_cs = dict()

for stepsize in [4e7, 5e7, 6e7, 7e7, 8e7, 9e7]:
    solution = tiht_mpa(A, y, 2 * RANK, stepsize=stepsize, batchsize=150)
    result_cs[stepsize] = list(a for a in Progress(RuntimeSlice(solution, 30), rettime=True))


    

    




100% |#####################################| Elapsed Time: 0:00:30 / 00:00:30
100% |#####################################| Elapsed Time: 0:00:30 / 00:00:30
100% |#####################################| Elapsed Time: 0:00:30 / 00:00:30
100% |#####################################| Elapsed Time: 0:00:30 / 00:00:30
100% |#####################################| Elapsed Time: 0:00:30 / 00:00:30
100% |#####################################| Elapsed Time: 0:00:30 / 00:00:30

In [11]:

    

pl.plot([t for t, _ in result], [mp.normdist(X, Xsharp) for _, (Xsharp, _) in result], label='adaptive')
for key, val in result_cs.items():
    pl.plot([t for t, _ in val], [mp.normdist(X, Xsharp) for _, Xsharp in val], label=key)

pl.ylim((0.95, 1.0))
pl.legend()


    

    
Out[11]:





<matplotlib.legend.Legend at 0x7fdf8d4da6a0>






    

In [15]:

    

res_mpa = [(mp.zero(SITES, LDIM, 1), None)]


    

In [16]:

    

ITSTEPS = 100
BATCHSIZE = 500


solution = tiht_mpa(A, y, 2 * RANK, stepsize='adaptive', retinfo=RETINFO,
                   batchsize=BATCHSIZE, X_init=res_mpa[-1][0])
with Timer():
    res_mpa += [p for p in islice(Progress(solution, max_value=ITSTEPS), ITSTEPS)]

pl.plot([mp.normdist(X, X_sharp) for X_sharp, _ in res_mpa], label='mpnum')


    

    




 28% ( 29 of 100) |######                  | Elapsed Time: 0:00:16 ETA: 0:00:36





    




Timer: 16.624500513076782s






    




 30% ( 30 of 100) |#######                 | Elapsed Time: 0:00:16 ETA: 0:00:35





    




---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-16-09152e03429a> in <module>()
      6                    batchsize=BATCHSIZE, X_init=res_mpa[-1][0])
      7 with Timer():
----> 8     res_mpa += [p for p in islice(Progress(solution, max_value=ITSTEPS), ITSTEPS)]
      9 
     10 pl.plot([mp.normdist(X, X_sharp) for X_sharp, _ in res_mpa], label='mpnum')

<ipython-input-16-09152e03429a> in <listcomp>(.0)
      6                    batchsize=BATCHSIZE, X_init=res_mpa[-1][0])
      7 with Timer():
----> 8     res_mpa += [p for p in islice(Progress(solution, max_value=ITSTEPS), ITSTEPS)]
      9 
     10 pl.plot([mp.normdist(X, X_sharp) for X_sharp, _ in res_mpa], label='mpnum')

/home/dsuess/Code/Pythonlibs/tools/helpers.py in __iter__(self)
    165         """Fetch next object from the iterable"""
    166         self.start()
--> 167         for n, val in enumerate(self._iterable):
    168             self.update(n)
    169             yield val

/files/home/part5/dsuess/Natty/Code/tiht.py in tiht_mpa(A, observation, bdim, X_init, stepsize, batchsize, retinfo, compargs)
    166         y_sharp = np.array([mpsp.inner_prod_mps(a, X) for a in A_])
    167         g = mp.sumup(A_, weights=y_ - y_sharp)
--> 168         g.compress(bdim=2 * bdim, **compargs)
    169 
    170         if stepsize == 'adaptive':

/files/home/part5/dsuess/Natty/Code/mpnum/mpnum/mparray.py in compress(self, method, **kwargs)
    800         """
    801         if method == 'svd':
--> 802             return self._compress_svd(**kwargs)
    803         elif method == 'swdsweep':
    804             pass

/files/home/part5/dsuess/Natty/Code/mpnum/mpnum/mparray.py in _compress_svd(self, bdim, relerr, direction)
    846 
    847         if direction == 'right':
--> 848             self.normalize(right=1)
    849             return self._compress_svd_r(bdim, relerr)
    850         elif direction == 'left':

/files/home/part5/dsuess/Natty/Code/mpnum/mpnum/mparray.py in normalize(self, left, right)
    671             self._lnormalize(lnormalize)
    672         if current_rnorm > rnormalize:
--> 673             self._rnormalize(rnormalize)
    674 
    675     def _lnormalize(self, to_site):

/files/home/part5/dsuess/Natty/Code/mpnum/mpnum/mparray.py in _rnormalize(self, to_site)
    712             # can be accounted by adapting bond dimension here
    713             self._ltens[site] = q.T.reshape((-1, ) + ltens.shape[1:])
--> 714             self._ltens[site - 1] = matdot(self._ltens[site - 1], r.T)
    715 
    716         self._lnormalized = min(to_site - 1, lnormal)

/files/home/part5/dsuess/Natty/Code/mpnum/mpnum/_tools.py in matdot(A, B, axes)
    105 def matdot(A, B, axes=((-1,), (0,))):
    106     """np.tensordot with sane defaults for matrix multiplication"""
--> 107     return np.tensordot(A, B, axes=axes)
    108 
    109 

/home/dsuess/local/conda/lib/python3.5/site-packages/numpy/core/numeric.py in tensordot(a, b, axes)
   1330     at = a.transpose(newaxes_a).reshape(newshape_a)
   1331     bt = b.transpose(newaxes_b).reshape(newshape_b)
-> 1332     res = dot(at, bt)
   1333     return res.reshape(olda + oldb)
   1334 

KeyboardInterrupt: 

In [ ]: